/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAWSWave
  \brief To calculate bound state radial wavefunction of a valence nucleon in
  a Woods-Saxon potential. Just for unit-test purposes.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022/03/07
  \date Last modified: 2022/03/07 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include "config.h"
#ifdef CATCH_CNOK
#include "catch.h"
#endif
#include <algorithm>
#include <iostream>
#include <cstring>
#include <fstream>
#include "TAWSWave.h"
#include "TAMath.h"
#include "TAInterpolate.h"
#include "TAIntegrate.h"

using std::max;
using std::ofstream;
using std::cout;
using std::endl;

#define polint TAInterpolateD::PolyInter
#define romb TAIntegrate<>::Romberg
#define sqr TAMath::sqr

// physical constants for calculation of the potential
static const double HBARC = TAMath::hbarc();
static const double RP2 = sqr(HBARC/TAMath::PionPMMass()); // pion compton wavelength^2
static const double e2 = HBARC*TAMath::FineStructureConstant(); // e^2/(4pi*epsilon0)

TAWSWave::TAWSWave() : TAODESolver(), fSolved(false), fxx(nullptr), fyy2(nullptr), fcnt(0){}

TAWSWave::~TAWSWave(){
  if(fxx) delete [] fxx;
  if(fyy2) delete [] fyy2;
} // end the dtor

void TAWSWave::derivs(double r, const double *y, double *dydx){
  dydx[0] = y[1];
  dydx[1] = (VT(r)-fE)*y[0];
} // end member function derivs

void TAWSWave::Bound(){
  if(fSolved) return;
  fSolved = true; // has been called

  SetEPS(1.e-6); // relative error at each step of RKstepper
  SetInitialStep(1.e-4);
  const double r0 = 1.e-10, rInf = 20.;
  SetMinStep(1.e-14);
  // the mindxsave is so that the x-range covers the whole interval //
  const int NMAX = 1000;
  SetSave(NMAX, (rInf-r0)/NMAX, 2); // (nmax, h1, nvar)
  n = 4; l = 0;
  // static const double EREF = -41.232607772;
  fE = -40.; // the initial guess for the eigenenergy
  double rm = 0.; // the turning point of total effective potential
  double dE = 1.e200;
  fxx = new double[NMAX*2]{}; // the concatenated solution: x
  fyy = new double[NMAX*2]{}; // the concatenated solution: y
  fyy2 = new double[NMAX*2]{}; // the concatenated solution: y^2
  // iteration to find the optimal eigenenergy //
  cout.precision(12); // DEBUG
  while(fabs(dE) > GetEPS()){ // usually E has the same accuracy with EPS
    ofstream fout("vis/rRl.txt");
    ofstream fout2("vis/rRl2.txt");
    // kinetic energy is postive at the inner part of the potential //
    if(fE < V(r0)) fE = V(r0) + fabs(V(r0)*0.01); // an emperical try
    // calculate the turning point rm as the first point from infinity where f(rm) = 0 //
    const double dr = rInf/NMAX;
    for(int i = NMAX; i--;) if(VT(dr*i)-fE < 0.){ rm = dr*i; break; }
    //-- solve the ODE --//
    double ri[NMAX]{}, ui[NMAX]{}; // the interior solution: r0->rm
    double re[NMAX]{}, ue[NMAX]{}; // the exterior solution: rInf->rm
    // integrate from r0 to rm //
    const double k = sqrt(fE-V(r0));
    const double x = r0*k, x2 = x*x, x4 = x2*x2, x6 = x4*x2;
    double l3 = (2.*l+3.), l5 = l3*(2.*l+5.), l7 = l5*(2.*l+7.);
    l3 = x2/(2.*l3); l5 = x4/(8.*l5); l7 = x6/(48.*l7);
    fy[0] = x*(1.-l3+l5-l7); // accuracy: O(x^8)*x
    fy[1] = k*(l+1.)-l3*(l+3.)+l5*(l+5.)-l7*(l+7.); // accuracy: O(x^8)
    ODEIntegrator(r0, rm);
    const int cnti = GetSolution(ri, ui); // extract the solution
    const double li = fy[1]/fy[0];
    // integrate from rInf to rm //
    const double beta = sqrt(fabs(fE));
    fy[0] = 1.7e-40; // exp(-beta*rInf)
    fy[1] = -beta*fy[0]; // -beta*
    ODEIntegrator(rInf, rm);
    const int cnte = GetSolution(re, ue);
    const double le = fy[1]/fy[0];
    // concatenate the solution //
    for(int i = 0; i < cnti; i++){
      fxx[i] = ri[i];
      fyy[i] = ui[i];
      fyy2[i] = sqr(fyy[i]);
    } // end for
    const double sca = ui[cnti-1]/ue[cnte-1];
    for(int i = 1; i < cnte; i++){
      const int i0 = cnti+i-1, i1 = cnte-i-1; // r=rm is dropped
      fxx[i0] = re[i1];
      fyy[i0] = ue[i1] * sca;
      fyy2[i0] = sqr(fyy[i0]);
    } // end for
    const int cnt = cnti+cnte-1; fcnt = cnt;
    const int nodes = TAMath::countNodes(fyy, cnt);
    // normalize the function //
    const double cc = pow(romb(cnt, fxx, fyy2), -0.5);
    for(int i = cnt; i--;){ fyy[i] *= cc; fyy2[i] = sqr(fyy[i]); }
    //////////// DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG //
    // cout << "li: " << li << " \tle: " << le << endl; // DEBUG
    // cout << "li-le: " << li-le << endl; // DEBUG
    for(int i = 0; i < cnt; i++) fout << fxx[i] << " " << fyy[i] << endl; // DEBUG
    for(int i = 0; i < cnt; i++) fout2 << fxx[i] << " " << fyy2[i] << endl; // DEBUG
    // cout << "Rl written to file" << endl; // DEBUG
    // cout << "nodes: " << nodes << endl; // DEBUG
    //////////// DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG //
    // cout << "______fE_before: " << fE << endl; // DEBUG
    if(nodes > n) fE *= 1.1;
    if(nodes < n) fE *= 0.9;
    if(nodes == n){ // calculate and apply the energy correction //
      dE = (li-le)*fyy2[cnti-1];
      fE += dE;
    } // end else
    // cout << "Eref: " << EREF << endl; // DEBUG
    // cout << "______fE_after: " << fE << endl; // DEBUG
    // cout << "dE: " << dE << endl; // DEBUG
    // cout << "|E-Eref|*1e9: " << fabs(fE-EREF)*1.e9 << endl; // DEBUG
    // getchar(); // DEBUG
  } // end while
  // cout << "\033[32;43;1m___BINGO___BINGO___BINGO\033[0m" << endl; exit(1); // DEBUG
  // getchar(); getchar(); getchar(); getchar(); // DEBUG
} // end member function Bound

/// \retval Note that Rl=u/r; u = Rl*r. This method returns Rl, i.e. u/r, NOT u
/// implement an interpolation
double TAWSWave::Rl(double r){ // u/r
  if(!fSolved) Bound();

  static const int NP = 4; // polynomial interpolation: Nof points
  if(r < 0.) TAException::Error("TAWSWave", "Rl: input r(%f) is minus.", r);
  if(r > fxx[fcnt-1]) TAException::Error("TAWSWave", "Rl: input r(%f) out of bounds.", r);
  return polint(fxx, fyy, fcnt, NP, r)/r;
} // end member function Rl



////////////////////////////// UNIT TEST ///////////////////////////////////////
// potential for the unit-test method Bound
double TAWSWave::V(double r) const{
  const double u0 = -50., a = 0.6, r0 = 7., u1 = -u0/a;
  const double t = exp((r-r0)/a);
  double v = u0/(1.+t) + u1*t/sqr(1.+t);
  return v; // *f2MuOverHbar2
} // end member function V
double TAWSWave::VT(double r) const{
  if(l){
    if(0. == r) TAException::Error("TAWSWave", "VT: r is 0 for nonzero l");
    return V(r) + l*(l+1.)/(r*r);
  } // end if(l)
  return V(r);
} // end member function VT

#ifdef CATCH_CNOK
TEST_CASE("Solving Radial equation with WS-Poential", "[ws]"){
  TAWSWave b;
  b.Bound();
  CHECK(b.E() == Approx(-41.232607772).epsilon(1.e-5));
} // end TEST_CASE ws
#endif
